	set more off
	global directory ""

	////////////////////////////////////////
	// 									  //
	// 				FIGURE 1	  		  //
	// 									  //
	////////////////////////////////////////
	
	global variables "compulsory_parent SLA_parent noqual_parent CSE_parent Olevel_parent Alevel_parent college_parent YoS_parent Inc_less_18k_parent Inc_more_31k_parent Inc_more_52k_parent Inc_more_100k_parent "
	global nvariables: word count $variables 
		
	use "$directory/Two_Generations_Sample", clear
		
	local i = 1
	foreach var in $variables {
		summ `var'
		replace `var' = (`var' - r(mean))/r(sd)
		clonevar depvar`i' = `var'
		clonevar score`i' = score
		local i = `i' + 1
	}
			
	forvalues i = 1/$nvariables {
			
		reg depvar`i' score`i', r
		estimates store parental`i'
				
		reg depvar`i' score`i' parentsscore, r
		estimates store siblings`i'

	}
		
	label variable score1	"Dropped at CSA" 
	label variable score2	"School-Leaving Age" 
	label variable score3	"No Qualifications" 
	label variable score4	"CSE"
	label variable score5	"O-Level"
	label variable score6	"A-Level" // 
	label variable score7	"College" 
	label variable score8	"Years of Schooling" 
	label variable score9	"Income < £18,000" 
	label variable score10	"Income > £31,000" 
	label variable score11	"Income > £52,000" 
	label variable score12	"Income > £100,000"
		

	coefplot (parental1 parental2 parental3 parental4 parental5 parental6 parental7 ///
			  parental8 parental9 parental10 parental11 parental12 ///
			 , mcolor(black) msymbol(O) ciopts(recast(rcap) color(black))) 					   /// 	
			 (siblings1 siblings2 siblings3 siblings4 siblings5 siblings6 siblings7  ///
			 siblings8 siblings9 siblings10 siblings11 siblings12 ///
			 , mcolor(red)   msymbol(T) ciopts(recast(rcap) color(red))) 						///
			 , keep(score1 score2 score3 score4 score5 score6 score7 score8 score9 score10 score11 score12) ///
			 headings(score1 = "{bf:Parental Education}" ///
			 score9 = "{bf:Parental Income}") ///
			 xline(0, lcolor(gs8) lwidth(thin)) scheme(s2mono) ///
			 graphregion(fcolor(white)) ylabel(, noticks) ///
			 xtitle("Effect of 1 SD Increase in Offspring's EA PGI ", size(medsmall)) /// 
			 coeflabels(, labsize(small)) ///
			 grid(glwidth(vthin)  within) /// 
			 legend(order(2 "Unconditional" 4 "Cond. on Parental PGI") position(3) ring(0) row(2) size(medsmall) region(lstyle(none))) //

	////////////////////////////////////////
	// 									  //
	// 				FIGURE 2	  		  //
	// 									  //
	////////////////////////////////////////
	
	use "$directory/Parents_Sample" if Inc_less_18k < ., clear
		
	twoway (kdensity residscore if Inc_less_18k == 1, lcolor(black) lpattern(dash)) (kdensity residscore if Inc_more_52k == 1, lcolor(red))   ///
		, ytitle("") title("Household Income", size(medsmall) ring(0) position(12)) /// 
		graphregion(fcolor(white)) ylabel(, noticks nogrid) xtitle("Residualized EA PGI") ///
		legend(order(1 "Household Income < £18,000" 2 "Household Income > £52,000") rows(1) region(lstyle(none))) 
		
		
	use "$directory/Parents_Sample" if home_ownership < ., clear
	
	twoway (kdensity residscore if home_ownership == 0, lcolor(black) lpattern(dash)) (kdensity residscore if home_ownership == 1, lcolor(red))   ///
		, ytitle("") title("Home Ownership", size(medsmall) ring(0) position(12)) /// 
		graphregion(fcolor(white)) ylabel(, noticks nogrid) xtitle("Residualized EA PGI") ///
		legend(order(1 "Does Not Own Living Accomodation" 2 "Owns Living Accomodation") rows(1) region(lstyle(none))) 
			 

	////////////////////////////////////////
	// 									  //
	// 				FIGURE 3	  		  //
	// 									  //
	////////////////////////////////////////
	
	use "$directory/Parents_Sample" if cond_wages < ., clear
		
	twoway (kdensity residscore if cond_wages <= 20000, lcolor(black) lpattern(dash)) (kdensity residscore if cond_wages >= 35000, lcolor(red))   ///
			, ytitle("") title("Occupational Wages", size(medsmall) ring(0) position(12)) /// 
			graphregion(fcolor(white)) ylabel(, noticks nogrid) xtitle("Residualized EA PGI") ///
			legend(order(1 "Occup. Wage < £20,000" 2 "Occup. Wage > £35,000") rows(1) region(lstyle(none))) 

	
	use "$directory/Parents_Sample" if onet_YoS_expert < ., clear
		
	
	twoway (kdensity residscore if onet_YoS_expert < 13, lcolor(black) lpattern(dash)) (kdensity residscore if onet_YoS_expert >= 16, lcolor(red))   ///
			, ytitle("") title("Occupation Educational Requirement", size(medsmall) ring(0) position(12)) /// 
			graphregion(fcolor(white)) ylabel(, noticks nogrid) xtitle("Residualized EA PGI") ///
			legend(order(1 "High School Degree or Less" 2 "Bachelor Degree or More") rows(1) region(lstyle(none))) 

		
	////////////////////////////////////////
	// 									  //
	// 				FIGURE 4	  		  //
	// 									  //
	////////////////////////////////////////
	
		
	use "$directory/Parents_Sample" if working == 1 & SOC2000 < ., clear
	
	replace SOC2000 = floor(SOC2000/10)
	
	foreach SOC2000 of numlist 111	112	113	114	115	116	117	118	121	122	123	211	212	213	221	231	232	241	242	243	244	245	311	312	313	321	322	323	331	341	342	343	344	351	352	353	354	355	356	411	412	413	414	415	421	511	521	522	523	524	531	532	541	542	543	549	611	612	613	621	622	623	629	711	712	721	811	812	813	814	821	822	911	912	913	914	921	922	923	924	925 {
		
		gen depvar = (SOC2000 == `SOC2000') 
			
		reg depvar score parentsscore, r
		
		gen coeff`SOC2000' = _b[score]
		
		test score
		gen pvalue`SOC2000' = r(p)
		
		drop depvar
	
	}
	
	gen n = _n
	keep if n == 1
	
	
	keep coeff* pvalue* 
	
	gen i = 1
	
	reshape long coeff pvalue, i(i) j(SOC2000)
	
	drop i
	
	
	merge 1:1 SOC2000 using "$directory/SOC2000_Labels.dta", nogen
	
	keep if SOC2000_label ~= ""
	
	
	sort coeff
	
	gen Z = 2 * _n
	
	scatter coeff Z, mlabel(SOC2000_label) msymbol(none) ///
	graphregion(fcolor(white)) ylabel(, nogrid)	///
	yline(0, lwidth(vthin) lcolor(gs8)) xtitle("Occupation") ///
	ytitle("Effect of EA PGI on Likelihood of Working as...") ///	
	xlabel(2(20)40, noticks labcolor(none)) xscale(range(2 40))  		
	
	////////////////////////////////////////
	// 									  //
	// 				FIGURE 5	  		  //
	// 									  //
	////////////////////////////////////////
	
	use "$directory/Parents_Sample" if compulsory < ., clear
		
	
	twoway (kdensity residscore if compulsory == 1, lcolor(black) lpattern(dash)) (kdensity residscore if compulsory == 0, lcolor(red))   ///
			, ytitle("") /// 
			graphregion(fcolor(white)) ylabel(, noticks nogrid) xtitle("Residualized EA PGI") ///
			legend(order(1 "Dropped at CSA" 2 "Stayed beyond CSA") rows(1) region(lstyle(none))) 
			
	////////////////////////////////////////
	// 									  //
	// 				FIGURE 6	  		  //
	// 									  //
	////////////////////////////////////////
	
	global educ "compulsory SLA noqual CSE Olevel Alevel college professional YoS NVQ"

	use "$directory/Parents_Sample", clear
	
	
	keep if lnwages < . & onet_YoS_expert < . & Inc_less_18k < . & YoS < . & SLA < .
	
	clonevar depvar1  = Inc_less_18k
	clonevar depvar2  = Inc_more_31k
	clonevar depvar3  = Inc_more_52k 
	clonevar depvar4  = Inc_more_100k 
	clonevar depvar5  = lnwages	
	clonevar depvar6  = onet_YoS_expert
	clonevar depvar7  = Inc_less_18k
	clonevar depvar8  = Inc_more_31k
	clonevar depvar9  = Inc_more_52k 
	clonevar depvar10 = Inc_more_100k 
	
	cap drop score2
	
	forvalues outcome = 1/10 {
		
		summ depvar`outcome'
		replace depvar`outcome' = (depvar`outcome' - r(mean))/r(sd)
		clonevar score`outcome'  = score

	}
		
	forvalues outcome = 1/4 {
			
		reg depvar`outcome' score`outcome' parentsscore if lnwages < ., r
		estimates store line`outcome'_uncond
				
		reg depvar`outcome' score`outcome' parentsscore lnwages, r 
		estimates store line`outcome'_cond
		
	}
	
		reg depvar5 score5 parentsscore if onet_YoS_expert < ., r
		estimates store line5_uncond
				
		reg depvar5 score5 parentsscore onet_YoS_expert, r 
		estimates store line5_cond

		
		reg depvar6 score6 parentsscore, r
		estimates store line6_uncond
				
		reg depvar6 score6 parentsscore $educ, r 
		estimates store line6_cond

		
	forvalues outcome = 7/10 {
			
		reg depvar`outcome' score`outcome' parentsscore, r
		estimates store line`outcome'_uncond
				
		reg depvar`outcome' score`outcome' parentsscore $educ, r 
		estimates store line`outcome'_cond
		
	}		
		
		label variable score1	"Control: Occupational Wages" 
		label variable score2	"Control: Occupational Wages" 
		label variable score3	"Control: Occupational Wages"
		label variable score4	"Control: Occupational Wages" 
		label variable score5	"Control: Occup. Educ. Requirement" 
		label variable score6	"Controls: Education" 
		label variable score7	"Controls: Education" 
		label variable score8	"Controls: Education" 
		label variable score9	"Controls: Education"
		label variable score10	"Controls: Education" 
			
		coefplot (line1_uncond, mcolor(black)  msymbol(O) ciopts(recast(rcap) color(black)) offset(0)) 	/// 	
				 (line1_cond, mcolor(red)    msymbol(T) ciopts(recast(rcap) color(red))) 	 /// 
				 (line2_uncond, mcolor(black)  msymbol(O) ciopts(recast(rcap) color(black))) 	/// 	
				 (line2_cond, mcolor(red)    msymbol(T) ciopts(recast(rcap) color(red)) ) 	 /// 
				 (line3_uncond, mcolor(black)  msymbol(O) ciopts(recast(rcap) color(black)) ) 	/// 	
				 (line3_cond, mcolor(red) 	  msymbol(T) ciopts(recast(rcap) color(red)) ) 	/// 	
				 (line4_uncond, mcolor(black)  msymbol(O) ciopts(recast(rcap) color(black)) ) 	 /// 
				 (line4_cond, mcolor(red)    msymbol(T) ciopts(recast(rcap) color(red)) ) 	 /// 
				 (line5_uncond, mcolor(black)  msymbol(O) ciopts(recast(rcap) color(black)) offset(0)) 	/// 	
				 (line5_cond, mcolor(red)   msymbol(T) ciopts(recast(rcap) color(red))) 	 /// 
				 (line6_uncond, mcolor(black)  msymbol(O) ciopts(recast(rcap) color(black))) 	/// 	
				 (line6_cond, mcolor(red)   msymbol(T) ciopts(recast(rcap) color(red)) ) 	 /// 
				 , keep(score1 score2 score3 score4 score5 score6 score7 score8 score9 score10) ///
				 graphregion(fcolor(white)) ylabel(, noticks) ///
				 grid(glwidth(vthin)  within) ///
				 coeflabels(, labsize(small)) xline(0, lcolor(gs8) lwidth(vthin)) ///
				xtitle("Effect of 1 SD Increase in the EA PGI") ///
				headings(score1 = "{bf:Dep. Var: Income < £18,000}" ///
				score2 = "{bf:Dep. Var: Income > £31,000}" ///
				score3 = "{bf:Dep. Var: Income > £52,000}" ///
				score4 = "{bf:Dep. Var: Income > £100,000}" ///
				score5 = "{bf:Dep. Var: Occupational Wages}" ///
				score6 = "{bf:Dep. Var: Occup. Educ. Requirement}" ) ///
				legend(order(2 "Without Controls" 20 "With Control(s)") row(1) size(medsmall) region(lstyle(none))) ///
				xscale(range(-0.04 0.08)) xlabel(-.04(.02)0.08)	
				
	////////////////////////////////////////
	// 									  //
	// 				FIGURE 7	  		  //
	// 									  //
	////////////////////////////////////////
	
	use "$directory/ROSLA_Sample" if residscore < ., clear
	
	
	drop residscore
	reg score parentsscore
	predict residscore, residuals
	
		cap drop temp
		gen temp = (year_birth*12 + month_birth) - (1957*12 + 9)
		gen YoB = floor(temp/12)
		
		bysort YoB: egen min = min(DoB)
		bysort YoB: egen max = max(DoB)
		gen x = (min + max)/2
		replace x = round(x)
	
	xtile pctile = residscore, nquantiles(2)
				
	reg CSE_Olevel c.DoB##c.DoB if DoB <  0 & pctile == 1 
	global pre_low_a = _b[_cons]
	global pre_low_b1 = _b[DoB]
	global pre_low_b2 = _b[c.DoB#c.DoB]

	reg CSE_Olevel c.DoB##c.DoB if DoB >= 0 & pctile == 1  
	global post_low_a = _b[_cons]
	global post_low_b1 = _b[DoB]
	global post_low_b2 = _b[c.DoB#c.DoB]

	reg CSE_Olevel c.DoB##c.DoB if DoB <  0 & pctile == 2
	global pre_high_a = _b[_cons]
	global pre_high_b1 = _b[DoB]
	global pre_high_b2 = _b[c.DoB#c.DoB]

	reg CSE_Olevel c.DoB##c.DoB if DoB >= 0 & pctile == 2 
	global post_high_a = _b[_cons]
	global post_high_b1 = _b[DoB]
	global post_high_b2 = _b[c.DoB#c.DoB]
				
	bysort YoB pctile: egen depvar = mean(CSE_Olevel)
				
	collapse (mean) CSE_Olevel (mean) x (mean) depvar, by(date_birth pctile) // (mean) YoB 
				
	gen DoB = date_birth - td(01sep1957)
							
	gen pre_low 	= $pre_low_a     + ($pre_low_b1 * DoB)    + ($pre_low_b2  * c.DoB#c.DoB) 	if DoB <= 0 & pctile == 1 
	gen post_low 	= $post_low_a    + ($post_low_b1 * DoB)   + ($post_low_b2  * c.DoB#c.DoB) 	if DoB >= 0 & pctile == 1				
	gen pre_high 	= $pre_high_a    + ($pre_high_b1 * DoB)   + ($pre_high_b2  * c.DoB#c.DoB) 	if DoB <= 0 & pctile == 2 
	gen post_high 	= $post_high_a   + ($post_high_b1 * DoB)  + ($post_high_b2  * c.DoB#c.DoB)	if DoB >= 0 & pctile == 2	


	bysort DoB pctile: gen n = _n
	replace x = . if n > 1
	
			
	twoway 	(scatter depvar    x if pctile == 1, msymbol(Oh)	mcolor(black) msize(vlarge)) ///
			(scatter depvar    x if pctile == 2, msymbol(X)	mcolor(red) msize(vlarge)) ///
			(line    pre_low   DoB if DoB <= 0 & pctile == 1, lcolor(black))  ///
			(line   post_low   DoB if DoB >= 0 & pctile == 1, lcolor(black))  ///
			(line    pre_high   DoB if DoB <= 0 & pctile == 2, lcolor(red))  ///
			(line   post_high   DoB if DoB >= 0 & pctile == 2, lcolor(red))  ///
			, ytitle("Fraction with an Educational Qualification") ///
			xtitle("Year of Birth") /// 
			graphregion(fcolor(white)) plotregion(lwidth(none)) ///
			scheme(s1mono) ///
			ylabel(, nogrid) ///
			title("`8'", size(medium) ring(0)) ///
			xlabel(-3470 "-10" -2740 "-8" -2009 "-6" -1279 "-4" -548 "-2" 182 "0" 913 "2" 1643 "4" 2374 "6"  3104 "8") ///
			legend(order(1 "Bottom Half" 2 "Top Half") rows(1)  region(lstyle(none))) //
						
	
	////////////////////////////////////////
	// 									  //
	// 				FIGURE 8	  		  //
	// 									  //
	////////////////////////////////////////

	use "$directory/Parents_Sample", clear
	
	global mediators "fluid neuro presbias"
	
	
	clonevar Y1 = SLA
	clonevar Y2 = noqual
	clonevar Y3 = college
	clonevar Y4 = YoS 

	
	forvalues outcome = 1/4 {
		summ Y`outcome'
		replace Y`outcome' = (Y`outcome' - r(mean))/r(sd)
	}
	
	
	ren fluid_intelligence fluid
	ren neuroticism neuro
	
		foreach mediator in $mediators {
			forvalues outcome = 1/4 {
				clonevar `mediator'_`outcome' = score
			}
		}

		
		forvalues outcome = 1/4 {				
				
			reg Y`outcome' fluid_`outcome' parentsscore if fluid < ., r
			estimates store fluid_`outcome'_1
				
			reg Y`outcome' fluid_`outcome' parentsscore fluid, r
			estimates store fluid_`outcome'_2
			
			
			reg Y`outcome' neuro_`outcome' parentsscore if neuro < ., r
			estimates store neuro_`outcome'_1
				
			reg Y`outcome' neuro_`outcome' parentsscore neuro, r
			estimates store neuro_`outcome'_2

			
			reg Y`outcome' presbias_`outcome' parentsscore if obese2 < . &  currently_smoking < ., r
			estimates store presbias_`outcome'_1
				
			reg Y`outcome' presbias_`outcome' parentsscore obese2 currently_smoking, r
			estimates store presbias_`outcome'_2

		}
		
	
		forvalues i = 1/4 {
			
			label variable neuro_`i'	"Control: Neuroticism" 
			label variable fluid_`i'	"Control: Fluid Intelligence" 
			label variable presbias_`i'	"Controls: Smoking & Obesity" 
			
		}
			
		coefplot (neuro_1_1, mcolor(black) msymbol(O) ciopts(recast(rcap) color(black)) offset(0)) 	/// 	
				 (neuro_1_2, mcolor(red)   msymbol(T) ciopts(recast(rcap) color(red))) 	 /// 
				 (fluid_1_1, mcolor(black) msymbol(O) ciopts(recast(rcap) color(black))) 	/// 	
				 (fluid_1_2, mcolor(red)   msymbol(T) ciopts(recast(rcap) color(red))) 	 /// 
				 (presbias_1_1, mcolor(black) msymbol(O) ciopts(recast(rcap) color(black))) 	/// 	
				 (presbias_1_2, mcolor(red)   msymbol(T) ciopts(recast(rcap) color(red))) 	 /// 
				 (neuro_2_1, mcolor(black) msymbol(O) ciopts(recast(rcap) color(black))) 	/// 	
				 (neuro_2_2, mcolor(red)   msymbol(T) ciopts(recast(rcap) color(red)) ) 	 /// 
				 (fluid_2_1, mcolor(black) msymbol(O) ciopts(recast(rcap) color(black))) 	/// 	
				 (fluid_2_2, mcolor(red)   msymbol(T) ciopts(recast(rcap) color(red))) 	 /// 
				 (presbias_2_1, mcolor(black) msymbol(O) ciopts(recast(rcap) color(black))) 	/// 	
				 (presbias_2_2, mcolor(red)   msymbol(T) ciopts(recast(rcap) color(red))) 	 /// 
				 (neuro_3_1, mcolor(black) msymbol(O) ciopts(recast(rcap) color(black)) ) 	/// 	
				 (neuro_3_2, mcolor(red) 	 msymbol(T) ciopts(recast(rcap) color(red)) ) 	/// 	
				 (fluid_3_1, mcolor(black) msymbol(O) ciopts(recast(rcap) color(black)) ) 	/// 	
				 (fluid_3_2, mcolor(red) msymbol(T) ciopts(recast(rcap) color(red)) ) 	/// 	
				 (presbias_3_1, mcolor(black) msymbol(O) ciopts(recast(rcap) color(black)) ) 	/// 	
				 (presbias_3_2, mcolor(red) msymbol(T) ciopts(recast(rcap) color(red)) ) 	/// 	
				 (neuro_4_1, mcolor(black) msymbol(O) ciopts(recast(rcap) color(black)) ) 	 /// 
				 (neuro_4_2, mcolor(red)   msymbol(T) ciopts(recast(rcap) color(red)) ) 	 /// 
				 (fluid_4_1, mcolor(black)   msymbol(O) ciopts(recast(rcap) color(black)) ) 	 /// 
				 (fluid_4_2, mcolor(red)   msymbol(T) ciopts(recast(rcap) color(red)) ) 	 /// 
				 (presbias_4_1, mcolor(black)   msymbol(O) ciopts(recast(rcap) color(black)) ) 	 /// 
				 (presbias_4_2, mcolor(red)   msymbol(T) ciopts(recast(rcap) color(red)) ) 	 /// 
				 , keep(fluid_* neuro_* presbias_*) ///
				 graphregion(fcolor(white)) ylabel(, noticks) ///
				 grid(glwidth(vthin)  within) ///
				 coeflabels(, labsize(small)) xline(0, lcolor(gs8) lwidth(vthin)) ///
				xtitle("Effect of 1 SD Increase in the EA PGI") ///
				headings(neuro_1 = "{bf:Dep. Var.: School-Leaving Age}" ///
				 neuro_2 = "{bf:Dep. Var.: No Qualification}" ///
				 neuro_3 = "{bf:Dep. Var.: College Degree}" ///
				 neuro_4 = "{bf:Dep. Var.: Years of Schooling}" ) ///
				legend(order(2 "Unconditional" 4 "Controlling for Potential Mediator") row(1) size(medsmall) region(lstyle(none))) //
				
	////////////////////////////////////////
	// 									  //
	// 				FIGURE 9	  		  //
	// 									  //
	////////////////////////////////////////
	
	use "$directory/Parents_Sample", clear
	
	keep if ts_33_0_0 >= td(01sep1957)
	
	twoway (kdensity residscore if started_smoking_age16 == 1, lcolor(black) lpattern(dash)) (kdensity residscore if started_smoking_age16 == 0, lcolor(red))   ///
			, ytitle("") /// 
			graphregion(fcolor(white)) ylabel(, noticks nogrid) xtitle("Residualized EA PGI") ///
			legend(order(1 "Started Smoking at Age 16 or Younger" 2 "Everyone Else") rows(1) region(lstyle(none))) 

	
	////////////////////////////////////////
	// 									  //
	// 			FIGURE 10A	  		      //
	// 									  //
	////////////////////////////////////////
	
	use "$directory/Parents_Sample", clear
	
	
	xtile terciles = parentsscore, nquantiles(3)
	
	gen low  = terciles == 1
	gen mid  = terciles == 2
	gen high = terciles == 3
	
	foreach tercile in low mid high {
		foreach var in score parentsscore {
			gen `var'_`tercile' = `var' * `tercile'
		}
	}
	
	cap drop x
	gen x = _n
	
	gen beta = .
	gen lb = .
	gen ub = .
	
	local x = 1
	foreach depvar in compulsory college { // noqual
		
		reg `depvar' score_high score_mid score_low parentsscore_high parentsscore_mid parentsscore_low high mid, r
		test (score_high = score_mid)
		test (score_high = score_low)
		test (score_low = score_mid)
		
		replace beta = _b[score_high] if x == `x' + 4
		replace beta = _b[score_mid]  if x == `x' + 2
		replace beta = _b[score_low]  if x == `x' 
	
		replace lb = _b[score_high] - (1.96 * _se[score_high]) if x == `x' + 4
		replace lb = _b[score_mid]  - (1.96 * _se[score_mid])  if x == `x' + 2
		replace lb = _b[score_low]  - (1.96 * _se[score_low])  if x == `x' 
	
		replace ub = _b[score_high] + (1.96 * _se[score_high]) if x == `x' + 4
		replace ub = _b[score_mid]  + (1.96 * _se[score_mid])  if x == `x' + 2
		replace ub = _b[score_low]  + (1.96 * _se[score_low])  if x == `x' 
		
		local x = `x' + 8
	
	}
	
	keep if x < 14
	
	twoway  (bar beta x if x == 1 | x == 9, color(ltblue)) ///
			(bar beta x if x == 3 | x == 11, color(emidblue)) ///
			(bar beta x if x == 5 | x == 13, color(edkblue)) ///
			(rcap ub lb x if x == 1 | x == 9,  lcolor(gs12)) ///
			(rcap ub lb x if x == 3 | x == 11, lcolor(gs13)) ///
			(rcap ub lb x if x == 5 | x == 13, lcolor(gs14)) ///
			, graphregion(fcolor(white)) plotregion(lwidth(none)) ///
			scheme(s1mono) ///
			legend(order (1 "Bottom Tercile" 2 "Middle Tercile" 3 "Top Tercile") rows(1)  region(lstyle(none))) ///
			text(.01  3 "...Dropping Out at CSA") ///
			text(.053 11 "...Graduating from College") ///
			xtitle("") ///
			xlabel("", nogrid noticks) ///
			yscale(range(-.05 .05)) ylabel(-0.05(.025)0.05, nogrid) ///
			ytitle("Effect of 1 SD Increase in EA PGI on Likelihood of...")
				
	////////////////////////////////////////
	// 									  //
	// 			FIGURE 10B	  		      //
	// 									  //
	////////////////////////////////////////

	use "$directory/Parents_Sample" if compulsory < . & college < ., clear
	
	cap drop x
	gen x = _n
	
	gen beta = .
	gen lb = .
	gen ub = .
	
	
	gen Y = .
		replace Y = 0	if compulsory == 1
		replace Y = 1	if compulsory == 0 & college == 0
		replace Y = 2	if college    == 1
	
	cap program drop boot_model
	program define boot_model, rclass
		
		oprobit Y score parentsscore, r
		
		gen Ylatent = (_b[score] * score) + (_b[parentsscore] * parentsscore) + rnormal()
		
		gen simulcompulsory = (Ylatent <=  _b[/cut1])
		gen simulcollege    = (Ylatent >   _b[/cut2])
			
		xtile terciles = parentsscore, nquantiles(3)
			
		gen low  = terciles == 1
		gen mid  = terciles == 2
		gen high = terciles == 3
		
		foreach tercile in low mid high {
			foreach var in score parentsscore {
				gen `var'_`tercile' = `var' * `tercile'
			}
		}		
		
		reg simulcompulsory score_high score_mid score_low parentsscore_high parentsscore_mid parentsscore_low high mid, r
		local bar1 = _b[score_low]
		local bar2 = _b[score_mid]
		local bar3 = _b[score_high]
		
		reg simulcollege score_high score_mid score_low parentsscore_high parentsscore_mid parentsscore_low high mid, r
		local bar4 = _b[score_low]
		local bar5 = _b[score_mid]
		local bar6 = _b[score_high]

		forvalues i = 1/6 {
			
			return scalar b`i' = `bar`i''
			
		}
		
		cap drop Ylatent simulcompulsory simulcollege terciles low mid high score_high score_mid score_low parentsscore_high parentsscore_mid parentsscore_low
		
	end
	
	bootstrap 	b1=r(b1) ///
				b2=r(b2) ///	
				b3=r(b3) ///
				b4=r(b4) ///	
				b5=r(b5) ///
				b6=r(b6) ///	
				, reps(500) seed(64835639) nodrop: boot_model
	
	local x = 1
	forvalues i = 1(3)4 {
		
		local j = `i' + 1
		local k = `i' + 2
		
		replace beta = _b[b`i']  if x == `x'
		replace beta = _b[b`j']  if x == `x' + 2
		replace beta = _b[b`k']  if x == `x' + 4 
	
		replace lb = _b[b`i']  - (1.96 * _se[b`i'])  if x == `x'
		replace lb = _b[b`j']  - (1.96 * _se[b`j'])  if x == `x' + 2
		replace lb = _b[b`k']  - (1.96 * _se[b`k'])  if x == `x' + 4 
	
		replace ub = _b[b`i']  + (1.96 * _se[b`i'])  if x == `x'
		replace ub = _b[b`j']  + (1.96 * _se[b`j'])  if x == `x' + 2
		replace ub = _b[b`k']  + (1.96 * _se[b`k'])  if x == `x' + 4 
		
		local x = `x' + 8
		
	}
	
	keep if x < 14
	
	twoway  (bar beta x if x == 1 | x == 9, color(ltblue)) ///
			(bar beta x if x == 3 | x == 11, color(emidblue)) ///
			(bar beta x if x == 5 | x == 13, color(edkblue)) ///
			(rcap ub lb x if x == 1 | x == 9,  lcolor(gs12)) ///
			(rcap ub lb x if x == 3 | x == 11, lcolor(gs13)) ///
			(rcap ub lb x if x == 5 | x == 13, lcolor(gs14)) ///
			, graphregion(fcolor(white)) plotregion(lwidth(none)) ///
			scheme(s1mono) ///
			legend(order (1 "Bottom Tercile" 2 "Middle Tercile" 3 "Top Tercile") rows(1)  region(lstyle(none))) ///
			text(.01  3 "...Dropping Out at CSA") ///
			text(.053 11 "...Graduating from College") ///
			xtitle("") ///
			xlabel("", nogrid noticks) ///
			yscale(range(-.05 .05)) ylabel(-0.05(.025)0.05, nogrid) //
